﻿#include <occmesh/api/filling.h>

FillingBuilder::FillingBuilder(const std::vector<std::tuple<TopoDS_Edge, TopoDS_Face>> &mapList)
{
    for (int i = 0; i < mapList.size(); i++)
    {
        double first, second;
        auto [e, f] = mapList[i];
        if (!f.IsNull())
        {
            m_boundaryData.boundarySurfaces.push_back(BRep_Tool::Surface(f));
            m_boundaryData.boundaryTypes.push_back(FillingType::G1);
        }
        else
        {
            m_boundaryData.boundarySurfaces.push_back(nullptr);
            m_boundaryData.boundaryTypes.push_back(FillingType::G0);
        }
        m_boundaryData.boundaryCurves.push_back(BRep_Tool::Curve(TopoDS::Edge(e), first, second));
        m_boundaryData.curveParamsRanges.push_back(std::make_tuple(first, second));
    }
}

TopoDS_Face FillingBuilder::Build(const FillingParamData &paramData)
{
    SampleAndParameterize();
    ConstructPCurves();

    int times = 0;
    double bestG0 = std::numeric_limits<double>::max();
    while (times++ < m_paramData.maxIter)
    {
        ConstructIntegralCurves();
        UpdateInitialSurface();
        ConstructOptimizedSystem();
        Reparameterize();
        ConstructPCurves();
        ModifyWeights();

        auto [G0, G1] = OCCUtils::MeasureDivation(m_resultData.surface, m_boundaryData.boundaryCurves,
                                                  m_boundaryData.boundarySurfaces);

        std::cout << "ave G0 error: " << G0 << std::endl;
        std::cout << "ave G1 error: " << G1 << std::endl;
        if (G0 < m_paramData.maxG0 || G0 < bestG0)
        {
            m_resultData.best = m_resultData.surface;
            bestG0 = std::min(G0, bestG0);

            if (bestG0 < m_paramData.maxG0 && G1 < m_paramData.maxG1)
                break;
        }
        // G0 误差增加则立即退出
        else
            break;
    }

    auto wire = OCCUtils::TrimmSurface(m_resultData.best, m_boundaryData.boundary2dCurves);
    if (wire.IsNull())
        return BRepBuilderAPI_MakeFace(m_resultData.best, 1e-16);
    else
        return BRepBuilderAPI_MakeFace(m_resultData.best, wire);
}

void FillingBuilder::SampleAndParameterize()
{
    int curveCount = m_boundaryData.boundaryCurves.size();
    const auto &boundaryCurves = m_boundaryData.boundaryCurves;
    m_boundaryData.samplePoints.clear();
    m_boundaryData.boundarySampleNums.clear();
    m_boundaryData.boundarySampleParams.clear();

    std::vector<double> curveLengths;
    for (int i = 0; i < curveCount; i++)
    {
        GeomAdaptor_Curve curveAdaptor(boundaryCurves[i]);
        double length = GCPnts_AbscissaPoint::Length(curveAdaptor, boundaryCurves[i]->FirstParameter(),
                                                     boundaryCurves[i]->LastParameter(), 1e-6);
        curveLengths.push_back(length);
    }
    double maxBoundaryCurveLength = *std::max_element(curveLengths.begin(), curveLengths.end());

    for (int i = 0; i < curveCount; i++)
    {
        int sampleNum =
            std::max(static_cast<int>(curveLengths[i] / maxBoundaryCurveLength * m_paramData.maxSampleNum), 3);
        m_boundaryData.boundarySampleNums.push_back(sampleNum);
    }

    for (int i = 0; i < curveCount; i++)
    {
        auto [tmin, tmax] = m_boundaryData.curveParamsRanges[i];
        auto sample = m_boundaryData.boundarySampleNums[i];
        auto dt = (tmax - tmin) / (sample - 1);
        for (int j = 0; j < sample; j++)
        {
            auto t = tmin + j * dt;
            m_boundaryData.boundarySampleParams.emplace_back(t);
            m_boundaryData.samplePoints.emplace_back(boundaryCurves[i]->Value(t));
        }
    }

    // 初始投影参数化
    m_boundaryData.uvParams = OCCUtils::ProjectParameterize(m_boundaryData.samplePoints);
    const auto &uvParams = m_boundaryData.uvParams;
    m_resultData.umax = std::max_element(uvParams.begin(), uvParams.end(), [](const auto &a, const auto &b) {
                            return a.X() < b.X();
                        })->X();
    m_resultData.vmax = std::max_element(uvParams.begin(), uvParams.end(), [](const auto &a, const auto &b) {
                            return a.Y() < b.Y();
                        })->Y();
    m_resultData.umin = std::min_element(uvParams.begin(), uvParams.end(), [](const auto &a, const auto &b) {
                            return a.X() < b.X();
                        })->X();
    m_resultData.vmin = std::min_element(uvParams.begin(), uvParams.end(), [](const auto &a, const auto &b) {
                            return a.Y() < b.Y();
                        })->Y();

    // 初始化节点向量
    m_resultData.uKnots.clear();
    m_resultData.vKnots.clear();

    double umax = m_resultData.umax;
    double umin = m_resultData.umin;
    double vmax = m_resultData.vmax;
    double vmin = m_resultData.vmin;
    for (int i = 0; i < m_paramData.uDegree + 1; i++)
        m_resultData.uKnots.push_back(umin - m_paramData.parameterFieldAdd * (umax - umin));
    for (int i = 0; i < m_paramData.uDegree + 1; i++)
        m_resultData.uKnots.push_back(umax + m_paramData.parameterFieldAdd * (umax - umin));

    for (int i = 0; i < m_paramData.vDegree + 1; i++)
        m_resultData.vKnots.push_back(vmin - m_paramData.parameterFieldAdd * (vmax - vmin));
    for (int i = 0; i < m_paramData.vDegree + 1; i++)
        m_resultData.vKnots.push_back(vmax + m_paramData.parameterFieldAdd * (vmax - vmin));
}

void FillingBuilder::ConstructPCurves()
{
    int sample = m_paramData.maxSampleNum;
    const auto &uvParams = m_boundaryData.uvParams;
    const auto &ts = m_boundaryData.boundarySampleParams;
    const auto &boundaryCurves = m_boundaryData.boundaryCurves;
    const auto &boundarySampleNums = m_boundaryData.boundarySampleNums;
    auto &boundary2dCurves = m_boundaryData.boundary2dCurves;
    boundary2dCurves.clear();

    int sampleId = 0;
    for (int i = 0; i < boundaryCurves.size(); i++)
    {
        std::vector<std::tuple<double, double, double>> tuv;
        for (int j = 0; j < boundarySampleNums[i]; j++)
        {
            tuv.emplace_back(ts[sampleId], uvParams[sampleId].X(), uvParams[sampleId].Y());
            sampleId++;
        }

        // 参数域上插值
        Handle(TColStd_HArray1OfReal) params = new TColStd_HArray1OfReal(1, tuv.size());
        Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d(1, tuv.size());
        for (int i = 0; i < tuv.size(); ++i)
        {
            auto [t, u, v] = tuv[i];
            params->SetValue(i + 1, t);
            points->SetValue(i + 1, gp_Pnt2d(u, v));
        }

        auto interp = Geom2dAPI_Interpolate(points, params, false, 1e-16);
        interp.Perform();
        boundary2dCurves.emplace_back(interp.Curve());
    }
}

void FillingBuilder::ConstructIntegralCurves()
{
    // 在积分节点处取值
    const auto &boundary2dCurves = m_boundaryData.boundary2dCurves;
    const auto &boundaryCurves = m_boundaryData.boundaryCurves;
    const auto &boundaryTypes = m_boundaryData.boundaryTypes;
    const auto &boundarySurfaces = m_boundaryData.boundarySurfaces;
    const auto &curveParamsRanges = m_boundaryData.curveParamsRanges;
    const auto &umin = m_resultData.umin, &umax = m_resultData.umax;
    const auto &vmin = m_resultData.vmin, &vmax = m_resultData.vmax;
    auto &integralNormals = m_integralData.integralNormals;
    auto &integralNodes = m_integralData.integralNodes;

    integralNormals.resize(boundaryCurves.size());
    integralNodes.resize(boundaryCurves.size());
    for (int i = 0; i < boundary2dCurves.size(); i++)
    {
        auto [tmin, tmax] = curveParamsRanges[i];
        double tlen = tmax - tmin;

        std::vector<std::tuple<double, double, double>> nodes;
        std::vector<gp_Vec> normals;
        for (int j = 0; j < OCCUtils::gaussNodes64.size(); j++)
        {
            double t = ((tmax - tmin) * OCCUtils::gaussNodes64[j] + tmax + tmin) / 2;
            gp_Pnt2d point2d = boundary2dCurves[i]->Value(t);

            // 更新参数范围
            // if (point2d.X() < umin)
            //     umin = point2d.X();
            // if (point2d.X() > umax)
            //     umax = point2d.X();
            // if (point2d.Y() < vmin)
            //     vmin = point2d.Y();
            // if (point2d.Y() > vmax)
            //     vmax = point2d.Y();
            nodes.emplace_back(t, point2d.X(), point2d.Y());

            if (boundaryTypes[i] == FillingType::G0)
                normals.emplace_back(0, 0, 0);
            else
            {
                // 节点法向
                gp_Vec normal;
                double u, v;
                auto point = boundaryCurves[i]->Value(t);
                OCCUtils::ProjectWithOCC(boundarySurfaces[i], point, u, v);

                gp_Vec D1U, D1V;
                boundarySurfaces[i]->D1(u, v, point, D1U, D1V);
                normal = D1U.Crossed(D1V);
                normal.Normalize();
                normals.emplace_back(normal);
            }
        }
        integralNormals[i] = std::move(normals);
        integralNodes[i] = std::move(nodes);
    }
}

void FillingBuilder::UpdateInitialSurface()
{
    const auto uDeg = m_paramData.uDegree, vDeg = m_paramData.vDegree;
    const auto &uKnots = m_resultData.uKnots, &vKnots = m_resultData.vKnots;
    auto &uN = m_paramData.uControlNum, &vN = m_paramData.vControlNum;

    // 根据生成的节点向量修改控制网格规模
    uN = uKnots.size() - uDeg - 2;
    vN = vKnots.size() - vDeg - 2;
    // 注意这里控制网格规模 -1 了

    // 构造初始曲面形状
    TColStd_Array1OfReal UKnotseq(1, uKnots.size());
    TColStd_Array1OfReal VKnotseq(1, vKnots.size());
    for (int i = 0; i < uKnots.size(); ++i)
        UKnotseq(i + 1) = uKnots[i];
    for (int i = 0; i < vKnots.size(); ++i)
        VKnotseq(i + 1) = vKnots[i];

    int uknotsLen = BSplCLib::KnotsLength(UKnotseq);
    int vknotsLen = BSplCLib::KnotsLength(VKnotseq);
    TColStd_Array1OfReal UKnots(1, uknotsLen), VKnots(1, vknotsLen);
    TColStd_Array1OfInteger UMults(1, uknotsLen), VMults(1, vknotsLen);
    BSplCLib::Knots(UKnotseq, UKnots, UMults);
    BSplCLib::Knots(VKnotseq, VKnots, VMults);

    TColgp_Array2OfPnt Poles(1, uN + 1, 1, vN + 1);
    m_resultData.surface = new Geom_BSplineSurface(Poles, UKnots, VKnots, UMults, VMults, uDeg, vDeg);
}

void FillingBuilder::ConstructOptimizedSystem()
{
    const auto uDeg = m_paramData.uDegree, vDeg = m_paramData.vDegree;
    const auto &uN = m_paramData.uControlNum, &vN = m_paramData.vControlNum;
    const auto &uKnots = m_resultData.uKnots, &vKnots = m_resultData.vKnots;
    const auto &boundary2dCurves = m_boundaryData.boundary2dCurves;
    const auto &boundaryCurves = m_boundaryData.boundaryCurves;
    const auto &boundaryTypes = m_boundaryData.boundaryTypes;
    const auto &boundarySurfaces = m_boundaryData.boundarySurfaces;
    const auto &curveParamsRanges = m_boundaryData.curveParamsRanges;
    const auto &integralNormals = m_integralData.integralNormals;
    const auto &integralNodes = m_integralData.integralNodes;
    const auto &weightG0 = m_paramData.weightG0, &weightG1 = m_paramData.weightG1;
    const auto &weightRocinBending = m_paramData.weightRocinBending;
    const auto &weightBending = m_paramData.weightBending;
    const auto &weightStretching = m_paramData.weightStretching;

    Eigen::MatrixXd A((uN + 1) * (vN + 1), (uN + 1) * (vN + 1));
    A.setZero();

    // 填充 A 矩阵
    std::vector<std::vector<double>> uders, vders;
    for (int iu = uDeg; iu <= uN; ++iu)
    {
        if ((uKnots[iu + 1] - uKnots[iu]) < 1e-12)
            continue;

        for (int iv = vDeg; iv <= vN; ++iv)
        {
            if ((vKnots[iv + 1] - vKnots[iv]) < 1e-12)
                continue;

            for (int i = 0; i <= 5; i++)
            {
                auto uders =
                    OCCUtils::DerivBasisFuncs(3, iu,
                                              OCCUtils::gaussNodes[4][i] * (uKnots[iu + 1] - uKnots[iu]) / 2.0 +
                                                  (uKnots[iu + 1] + uKnots[iu]) / 2.0,
                                              uDeg, uKnots);
                for (int j = 0; j <= 5; j++)
                {
                    double coeff = OCCUtils::gaussWeights[4][i] * OCCUtils::gaussWeights[4][j] *
                                   (uKnots[iu + 1] - uKnots[iu]) * (vKnots[iv + 1] - vKnots[iv]) / 4.0;
                    auto vders =
                        OCCUtils::DerivBasisFuncs(3, iv,
                                                  OCCUtils::gaussNodes[4][j] * (vKnots[iv + 1] - vKnots[iv]) / 2.0 +
                                                      (vKnots[iv + 1] + vKnots[iv]) / 2.0,
                                                  vDeg, vKnots);
                    for (int i1 = 0; i1 <= uDeg; ++i1)
                    {
                        for (int i2 = 0; i2 <= uDeg; ++i2)
                        {
                            for (int l1 = 0; l1 <= vDeg; ++l1)
                            {
                                for (int l2 = 0; l2 <= vDeg; ++l2)
                                {
                                    A((iu - uDeg + i1) * (vN + 1) + (iv - vDeg + l1),
                                      (iu - uDeg + i2) * (vN + 1) + (iv - vDeg + l2)) +=
                                        ((uders[3][i1] * uders[3][i2] * vders[0][l1] * vders[0][l2] +
                                          uders[2][i1] * uders[2][i2] * vders[1][l1] * vders[1][l2] * 3 +
                                          uders[1][i1] * uders[1][i2] * vders[2][l1] * vders[2][l2] * 3 +
                                          uders[0][i1] * uders[0][i2] * vders[3][l1] * vders[3][l2]) *
                                             weightRocinBending +
                                         (uders[2][i1] * uders[2][i2] * vders[0][l1] * vders[0][l2] +
                                          uders[1][i1] * uders[1][i2] * vders[1][l1] * vders[1][l2] * 2 +
                                          uders[0][i1] * uders[0][i2] * vders[2][l1] * vders[2][l2]) *
                                             weightBending +
                                         (uders[1][i1] * uders[1][i2] * vders[0][l1] * vders[0][l2] +
                                          uders[0][i1] * uders[0][i2] * vders[1][l1] * vders[1][l2]) *
                                             weightStretching) *
                                        coeff;
                                }
                            }
                        }
                    }
                }
            }
        }
    }

    Eigen::MatrixXd D0((uN + 1) * (vN + 1), (uN + 1) * (vN + 1));
    D0.setZero();
    Eigen::VectorXd b0x((uN + 1) * (vN + 1));
    b0x.setZero();
    Eigen::VectorXd b0y((uN + 1) * (vN + 1));
    b0y.setZero();
    Eigen::VectorXd b0z((uN + 1) * (vN + 1));
    b0z.setZero();

    // 填充 D0, b0x, b0y, b0z 矩阵
    for (int i = 0; i < boundaryCurves.size(); ++i)
    {
        auto &tuvs = integralNodes[i];
        auto &curve = boundaryCurves[i];
        auto [tmin, tmax] = curveParamsRanges[i];
        double dt = (tmax - tmin) / 2;
        for (int j = 0; j < tuvs.size(); ++j)
        {
            auto [t, u, v] = tuvs[j];
            int uspan = OCCUtils::FindSpan(u, uDeg, uKnots);
            int vspan = OCCUtils::FindSpan(v, vDeg, vKnots);

            auto uval = OCCUtils::BasisFuncs(uspan, u, uDeg, uKnots);
            auto vval = OCCUtils::BasisFuncs(vspan, v, vDeg, vKnots);

            auto point = curve->Value(t);
            // calculate position deviation and integrate
            for (int i1 = 0; i1 <= uDeg; ++i1)
            {
                for (int l1 = 0; l1 <= vDeg; ++l1)
                {
                    auto temp = uval[i1] * vval[l1] * OCCUtils::gaussWeights64[j] * dt;
                    b0x((uspan - uDeg + i1) * (vN + 1) + (vspan - vDeg + l1)) += point.X() * temp;
                    b0y((uspan - uDeg + i1) * (vN + 1) + (vspan - vDeg + l1)) += point.Y() * temp;
                    b0z((uspan - uDeg + i1) * (vN + 1) + (vspan - vDeg + l1)) += point.Z() * temp;

                    for (int i2 = 0; i2 <= uDeg; ++i2)
                    {
                        for (int l2 = 0; l2 <= vDeg; ++l2)
                        {
                            D0((uspan - uDeg + i1) * (vN + 1) + (vspan - vDeg + l1),
                               (uspan - uDeg + i2) * (vN + 1) + (vspan - vDeg + l2)) += (uval[i2] * vval[l2]) * temp;
                        }
                    }
                }
            }
        }
    }

    Eigen::MatrixXd D1xx((uN + 1) * (vN + 1), (uN + 1) * (vN + 1));
    D1xx.setZero();
    Eigen::MatrixXd D1xy((uN + 1) * (vN + 1), (uN + 1) * (vN + 1));
    D1xy.setZero();
    Eigen::MatrixXd D1xz((uN + 1) * (vN + 1), (uN + 1) * (vN + 1));
    D1xz.setZero();
    Eigen::MatrixXd D1yy((uN + 1) * (vN + 1), (uN + 1) * (vN + 1));
    D1yy.setZero();
    Eigen::MatrixXd D1yz((uN + 1) * (vN + 1), (uN + 1) * (vN + 1));
    D1yz.setZero();
    Eigen::MatrixXd D1zz((uN + 1) * (vN + 1), (uN + 1) * (vN + 1));
    D1zz.setZero();

    // 填充 D1xx, D1xy, D1xz, D1yy, D1yz, D1zz 矩阵
    for (int i = 0; i < boundaryCurves.size(); ++i)
    {
        if (boundaryTypes[i] == FillingType::G0)
            continue;

        auto &tuvs = integralNodes[i];
        auto &curve2d = boundary2dCurves[i];
        auto [tmin, tmax] = curveParamsRanges[i];
        double dt = (tmax - tmin) / 2;
        for (int j = 0; j < tuvs.size(); ++j)
        {
            auto [t, u, v] = tuvs[j];
            int uspan = OCCUtils::FindSpan(u, uDeg, uKnots);
            int vspan = OCCUtils::FindSpan(v, vDeg, vKnots);

            auto uval = OCCUtils::DerivBasisFuncs(2, uspan, u, uDeg, uKnots);
            auto vval = OCCUtils::DerivBasisFuncs(2, vspan, v, vDeg, vKnots);

            double ut, vt;
            gp_Vec2d tangent;
            gp_Pnt2d temp;
            curve2d->D1(t, temp, tangent);
            tangent.Normalize();
            ut = tangent.X();
            vt = tangent.Y();

            gp_Vec vec = integralNormals[i][j];
            // calculate normal deviation and integrate
            for (int i1 = 0; i1 <= uDeg; ++i1)
            {
                for (int i2 = 0; i2 <= uDeg; ++i2)
                {
                    for (int l1 = 0; l1 <= vDeg; ++l1)
                    {
                        for (int l2 = 0; l2 <= vDeg; ++l2)
                        {
                            double temp = (vt * uval[1][i1] * vval[0][l1] - ut * uval[0][i1] * vval[1][l1]) *
                                          (vt * uval[1][i2] * vval[0][l2] - ut * uval[0][i2] * vval[1][l2]) * dt *
                                          OCCUtils::gaussWeights64[j];

                            D1xx((uspan - uDeg + i1) * (vN + 1) + (vspan - vDeg + l1),
                                 (uspan - uDeg + i2) * (vN + 1) + (vspan - vDeg + l2)) += vec.X() * vec.X() * temp;
                            D1xy((uspan - uDeg + i1) * (vN + 1) + (vspan - vDeg + l1),
                                 (uspan - uDeg + i2) * (vN + 1) + (vspan - vDeg + l2)) += vec.X() * vec.Y() * temp;
                            D1xz((uspan - uDeg + i1) * (vN + 1) + (vspan - vDeg + l1),
                                 (uspan - uDeg + i2) * (vN + 1) + (vspan - vDeg + l2)) += vec.X() * vec.Z() * temp;
                            D1yy((uspan - uDeg + i1) * (vN + 1) + (vspan - vDeg + l1),
                                 (uspan - uDeg + i2) * (vN + 1) + (vspan - vDeg + l2)) += vec.Y() * vec.Y() * temp;
                            D1yz((uspan - uDeg + i1) * (vN + 1) + (vspan - vDeg + l1),
                                 (uspan - uDeg + i2) * (vN + 1) + (vspan - vDeg + l2)) += vec.Y() * vec.Z() * temp;
                            D1zz((uspan - uDeg + i1) * (vN + 1) + (vspan - vDeg + l1),
                                 (uspan - uDeg + i2) * (vN + 1) + (vspan - vDeg + l2)) += vec.Z() * vec.Z() * temp;
                        }
                    }
                }
            }
        }
    }

    Eigen::MatrixXd C((uN + 1) * (vN + 1) * 3, (uN + 1) * (vN + 1) * 3);

    C.block(0, 0, (uN + 1) * (vN + 1), (uN + 1) * (vN + 1)) = A + D0 * weightG0 + D1xx * weightG1;
    C.block(0, (uN + 1) * (vN + 1), (uN + 1) * (vN + 1), (uN + 1) * (vN + 1)) = D1xy * weightG1;
    C.block(0, 2 * (uN + 1) * (vN + 1), (uN + 1) * (vN + 1), (uN + 1) * (vN + 1)) = D1xz * weightG1;
    C.block((uN + 1) * (vN + 1), 0, (uN + 1) * (vN + 1), (uN + 1) * (vN + 1)) = D1xy * weightG1;
    C.block((uN + 1) * (vN + 1), (uN + 1) * (vN + 1), (uN + 1) * (vN + 1), (uN + 1) * (vN + 1)) =
        A + D0 * weightG0 + D1yy * weightG1;
    C.block((uN + 1) * (vN + 1), 2 * (uN + 1) * (vN + 1), (uN + 1) * (vN + 1), (uN + 1) * (vN + 1)) = D1yz * weightG1;
    C.block(2 * (uN + 1) * (vN + 1), 0, (uN + 1) * (vN + 1), (uN + 1) * (vN + 1)) = D1xz * weightG1;
    C.block(2 * (uN + 1) * (vN + 1), (uN + 1) * (vN + 1), (uN + 1) * (vN + 1), (uN + 1) * (vN + 1)) = D1yz * weightG1;
    C.block(2 * (uN + 1) * (vN + 1), 2 * (uN + 1) * (vN + 1), (uN + 1) * (vN + 1), (uN + 1) * (vN + 1)) =
        A + D0 * weightG0 + D1zz * weightG1;

    Eigen::VectorXd b(3 * (uN + 1) * (vN + 1));
    b.block(0, 0, (uN + 1) * (vN + 1), 1) = b0x * weightG0;
    b.block((uN + 1) * (vN + 1), 0, (uN + 1) * (vN + 1), 1) = b0y * weightG0;
    b.block(2 * (uN + 1) * (vN + 1), 0, (uN + 1) * (vN + 1), 1) = b0z * weightG0;

    // 求解系统
    Eigen::VectorXd x(3 * (uN + 1) * (vN + 1));
    x = C.llt().solve(b);

    auto const &px = x.block(0, 0, (uN + 1) * (vN + 1), 1);
    auto const &py = x.block((uN + 1) * (vN + 1), 0, (uN + 1) * (vN + 1), 1);
    auto const &pz = x.block(2 * (uN + 1) * (vN + 1), 0, (uN + 1) * (vN + 1), 1);

    // 构造结果
    for (int i = 0; i < uN + 1; i++)
    {
        for (int j = 0; j < vN + 1; j++)
        {
            gp_Pnt point(px(i * (vN + 1) + j, 0), py(i * (vN + 1) + j, 0), pz(i * (vN + 1) + j, 0));
            m_resultData.surface->SetPole(i + 1, j + 1, point);
        }
    }

    std::cout << "G0 Weight: " << weightG0 << " G1 Weight: " << weightG1 << std::endl;
    std::cout << "weightRocinBending: " << weightRocinBending << " weightBending: " << weightBending
              << " weightStretching: " << weightStretching << std::endl;
    std::cout << "Surface Degree: " << uDeg << " " << vDeg << std::endl;
    std::cout << "Surface Knots: " << m_resultData.surface->NbUKnots() << " " << m_resultData.surface->NbVKnots()
              << std::endl;
    std::cout << "Surface Poles: " << m_resultData.surface->NbUPoles() << " " << m_resultData.surface->NbVPoles()
              << std::endl;
}

void FillingBuilder::Reparameterize()
{
    auto &umax = m_resultData.umax, &umin = m_resultData.umin;
    auto &vmax = m_resultData.vmax, &vmin = m_resultData.vmin;

    // u, v 插入点数
    double ratio = (umax - umin) / (vmax - vmin);
    int uInsertNum = std::min(2.0, ratio <= 1 ? 1 : ratio);
    int vInsertNum = std::min(2.0, ratio >= 1 ? 1 : 1 / ratio);

    const auto &points = m_boundaryData.samplePoints;
    auto &uvParams = m_boundaryData.uvParams;
    std::vector<std::tuple<double, double, double>> uvd;
    for (int i = 0; i < points.size(); i++)
    {
        double u = uvParams[i].X(), v = uvParams[i].Y();
        double distance = OCCUtils::ProjectWithNearest(m_resultData.surface, points[i], u, v, 1e-6);
        uvParams[i] = gp_Pnt2d(u, v);
        uvd.push_back({u, v, distance});
    }

    std::sort(uvd.begin(), uvd.end(), [](const auto &a, const auto &b) { return std::get<2>(a) > std::get<2>(b); });

    auto &uKnots = m_resultData.uKnots;
    auto &vKnots = m_resultData.vKnots;

    const auto uDeg = m_paramData.uDegree, vDeg = m_paramData.vDegree;
    while (uInsertNum-- > 0)
    {
        auto uInsertPosition = OCCUtils::FindSpan(std::get<0>(uvd[uInsertNum]), uDeg, uKnots);
        auto uInsert = (uKnots.at(uInsertPosition) + uKnots.at(uInsertPosition + 1)) / 2;
        for (int i = 0; i < m_paramData.insertNum; i++)
            uKnots.insert(uKnots.begin() + uInsertPosition + 1, uInsert);
    }

    while (vInsertNum-- > 0)
    {
        auto vInsertPosition = OCCUtils::FindSpan(std::get<1>(uvd[vInsertNum]), vDeg, vKnots);
        auto vInsert = (vKnots.at(vInsertPosition) + vKnots.at(vInsertPosition + 1)) / 2;
        for (int i = 0; i < m_paramData.insertNum; i++)
            vKnots.insert(vKnots.begin() + vInsertPosition + 1, vInsert);
    }
}

void FillingBuilder::ModifyWeights()
{
    // m_paramData.weightBending /= 10;
    // m_paramData.weightStretching /= 10;
    m_paramData.weightG0 *= 40;
    m_paramData.weightG1 *= 20;
}
